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SED FITTING WITH MARKOV CHAIN MONTE CARLO: 
METHODOLOGY AND APPLICATION TO Z=3.1 LYa-EMITTING GALAXIES 
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ABSTRACT 

We present GalMC, a MCMC algorithm designed to fit the spectral energy distributions (SED) of 
galaxies to infer physical properties such as age, stellar mass, dust reddening, metallicity, redshift, and 
star formation rate. We describe the features of the code and the extensive tests conducted to ensure 
that our procedure leads to unbiased parameter estimation and accurate evaluation of uncertainties. 
We compare its performance to grid-based algorithms, showing that the efficiency in CPU time is ~ 
100 times better for MCMC for a three dimensional parameter space and increasing with the number 
of dimensions. We use GalMC to fit the stacked SEDs of two samples of Lyman Alpha Emitters 
(LAEs) at redshift z=3.1. Our fit reveals that the typical LAE detected in the IRAC 3.6 /im band 
has age = 0.67 [0.37 -1.81] Gyr and stellar mass = 3.2 [2.5 - 4.2] xlO^ Mq, while the typical LAE 
not detected at 3.6 ^m has age = 0.06 [0.01-0.2] Gyr and stellar mass = 2 [1.1 - 3.4] xlO® Mq. The 
SEDs of both stacks are consistent with the absence of dust. The data do not significantly prefer 
exponential with respect to constant star formation history. The stellar populations of these two 
samples are consistent with the previous study by Lai et al, with some differences due to the improved 
modeling of the stellar populations. A constraint on the metallicity of z=3.1 LAEs from broad-band 
photometry, requiring Z < Zq at 95% confidence, is found here for the first time. 



1. INTRODUCTION 

A galaxy's Spectral Energy Distribution (SED) con- 
tains information about its stellar population age, mass, 
star formation rate (SFR), dust content, and metallicity. 
Different physical processes leave their imprint in dif- 
ferent parts of the spectrum; the wider the wavelength 
coverage, the more robust the interpretation of the var- 
ious features of the SED is in terms of galaxy proper- 
ties. The rest-frame ultraviolet (UV), optical, and in- 
frared (IR) parts of the spectrum offer the combination 
of depth and angular resolution needed to obtain individ- 
ual photometry for all galaxies detected in a particular 
band, which is not generally possible in the far-infrared 
through sub-millimeter wavelengths where re-emission 
by dust dominates the luminosity. Spectroscopic infor- 
mation is generally required to determine metallicity and 
further probes of chemical enrichment, but spectroscopy 
of unbiased samples of dim galaxies is very difficult to 
obtain. This leaves photometric UV-through-IR SEDs 
as the most readily available probe of galaxy properties. 

SED fitting is the procedure of compa ring models to 
the observed galaxy SED (for reviews, see IGawised[200l 
iWalcher et al.ll2010l and references therein). Since the 
physical properties of the models are known, those of 
the data can be derived by maximizing the resemblance 
between data and models. The success and reliabil- 
ity of this method depend on the quality of the avail- 
able template spectra and on the robustness of the fit- 
ting algorithm. The ingredients of the spectral tem- 
plates typically are libraries of stellar spectra and sets 
of evolutionary tracks, allowing one to compute the ini- 
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tial spectrum of a collection of stars of different mass 
and to follow its evolution with time. Measuring and 
calibrating these quantities is not an easy process, and 
the available models cannot uniformly cover stellar pop- 
ulations of all ages, masses, or chemical content. The 
large variance intrinsic in the spectral templates is re- 
flected in the large number of available stellar popula- 
tion synthesis (SPS) codes, including but no t limited 
to PEGASE ((Fioc & Rocca-Volmcrangc"'1997), Maras- 
1998, 2005), GRASIL (Sil va et al. 1998.1. 
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Starb urst99 ( Leitherer et al.iri999t IVazquez fc Leithered 
l2005h . GALAXEV (Bruzual and Chario t 2003, Char- 
lot a nd Bruzual 2010), and GALEV ()Kotulla et al.l 
l2009f ). The scatter among the predictions of different 
SPS models is in itself an important source of sy stem- 
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atic uncert ainty (see e.g., iKannappan fc Gawiser 
iConrov fc~ G unn 2 010l and our discussion in Sec. |4]). 

Often, the comparison between model and data is 
made through a minimization, which provides us 
with a best-fit model. If the probability distribu- 
tion of th e parameters is close to a Gaussian {e.g., 
iLarson et~a l. 2010), the best fit is a good estimate 
of the expectation value for each parameter, and the 
corresponding uncertainties can be evaluated from 
relative likelihoods assuming a Gaussian profile. This 
approach can, however, be very dangerous in astron- 
omy. In fact, because of the unprecedented volume of 
available data, we can now hope to explore new aspects 
of the physics of galaxies, and this exciting perspective 
comes at the price of not having prior knowledge of 
the probability distribution we set out to investigate. 
Furthermore, the high degree of correlation between 
astrophysical processes makes probability distributions 
highly non-Gaussian. In a multi-dimensional space with 
heavy parameter degeneracies, the "best fit" spectrum 
is an overly aggressive compression of the information 
available in the SED, and the assumption of a Gaussian 
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likelihood is unjustified. As a result, it is necessary 
to switch to more sophisticated and robust statistical 
analysis techniques, such as Markov Chain Monte 
Carlo (MCMC) algorithms. While the use of simple 
minimization techniques is still very common, MCMC is 
becoming a popular sta tistical method within the SEP 
fittin g comniunity (e.g.. iPanter et al.ll2003t iSaiina et al.l 
20061 INilsson et all [20071 I201tt IConrov fc GunnI 120101 



Serra et al.ll2011l: iPirzkal et al.ll201lL Johnson et al 2011 



in prep); however, none of these algorithms is publicly 
available yet. Our MCMC code, GalMC, is available at 
|http://www.physics.rutgers.edu/~vacquaviva/web/GalMC.html| 

This paper is organized as follows. In Sec. [2]we present 
the general features of MCMC algorithms and describe in 
the detail GalMC, our MCMC code for SED fitting. Sec. 
[3] illustrates the extensive tests we performed on Mock 
catalogs to ensure that the algorithm works correctly and 
that our inferred parameter values and uncertainties are 
trustworthy. In Sec. |4]we present the results we obtain 
for two samples of Lyman Alpha Emitters (LAE) at red- 
shift z = 3.1, for various assumptions of star formation 
histories (SFH) and stellar populations synthesis (SPS) 
models; we also compare our results to those obtained for 
the same samples, but with a different technique, by Lai 
et al 2008 (hereafter LOS). We summarize our findings 
in Sec. [5] 

2. THE MCMC ALGORITHM FOR SED FITTING 

The process of parameter fitting involves two essential 
steps. First, we want to know what is the most repre- 
sentative value, or estimate, of the true value of each pa- 
rameter included in our fit. Second, we want to evaluate 
the uncertainties, i.e., how likely it is that the true value 
lies in a given interval in the vicinity of that estimate. 
In order to accomplish this task, one needs to determine 
the probability density function of the parameters, p(x), 
which can be used to compute the expectation value of 
any function of interest as 

(/(x)) = j dxi...da;„/(x)p(x). (1) 

In many cases where the probability density function is 
complicated and the dimensionality n of the parameter 
space is high, finding alternate ways of evaluating the 
above integral is required. This can achieved by informed 
sampling of the parameter space. A sample is a vector of 
coordinates in parameter space; the objective of the sam- 
pling process is to obtain a set of R independent samples 
Ti whose distribution is identical to that of the prob- 
ability density function p(x). In practice, this can be 
achieved if the density of sampled points is proportional 
to p(x). Once this collection of samples is obtained, it 
will by definition satisfy the property 
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for any function /, and performing integrals of the form 
of Eq. ([Ij becomes a trivial matter. Often, we will be 
interested in recovering the posterior probability distri- 
bution of parameters. According to Bayesian statistics, 
which we adopt, this depends on the product of the likeli- 
hood of observing the data, if that model were the truth. 



and on our beliefs on the distribution of parameters, or 
priors. 

A flexible method allowing one to create samples with 
minimal assumptions on the underlying probability dis- 
tribution of the parameters with respect to th e data is 
to us e Monte Carlo M arkov Chains (see e. g.. iMacKav 
2001 iNeail 119 93: Lew is fc Bridld I20q2at IVerde et al. 
2003: lHaiianl " l2007; .Hobson et al.ll2010[) . Monte Carlo 



methods are computational techniques that make use of 
random numbers. In this context, a random number is 
involved in the process of moving from one point in pa- 
rameter space to the next one. This succession of steps 
is called a chain. Because the transition to the next step 
only depends on the current point, having no memory of 
the previous points, the chain is called a Markov chain. 

How can we obtain the desired independent samples 
from a Markov chain? The chain needs to be built satis- 
fying certain rules that restrict the transition probability 
r(x, y), the probability of moving from a point x to t he 
point y. Two things are required {e.g., lMacKavll2003[ ): 

1. The desired probability distribution p is an invari- 
ant distribution of the chain: 



P(x) ^ J d'Y r(x, y) p{y). 



(3) 



2. The chain is ergodic, i.e., the probability distribu- 
tion of the chain tends to an invariant distribution 
of the chain, no matter what the initial conditions 
are. 

If these conditions are satisfied, after an initial period 
where the walk of the chain depends on the initial condi- 
tions (the burn-in phase), the chain will start sampling 
from the desired probability distribution (the posterior) . 
However, subsequent steps of the chain are strongly cor- 
related with each other. It is necessary to select points 
which are distant enough from each other (thin the chain) 
in order to achieve the independent samples needed in 
Eq. H 

There are two main tasks associated with the creation 
of chain. The first is how to explore the parameter space, 
namely, how to go efficiently from one point to the next 
one and how to decide whether or not the step just taken 
will become part of the chain. The latter process gen- 
erally involves evaluating how similar the model corre- 
sponding to a point in parameter space is to the data 
one wants to fit, based on the model's likelihood given 
the observed data, and on the chosen priors. The sec- 
ond task is to compute the likelihood of a given model at 
each point visited by the chain; we do so by using stel- 
lar population synthesis models to predict the observable 
quantities as a function of the input parameters. These 
two processes are described in detail below. 

2.1. Description and features of sampling algorithm 

The structure of our code, GalMC, is deeply in- 
debted to the publicly available algorithm CosmoMC 
(Lewis & Bridle 2002a). We use a Metropolis sampling 
algorithm (Hastings 1970), where a trial step Xq -|- Ax 
from the current position in parameter space Xq is ac- 
cepted with probability 
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where p(xo) is the posterior probability evaluated at xq. 

The trial step Xq + Ax is a random deviate drawn 
from a distribution, (7(x), known as the proposal density. 
In principle, Metropolis sampling converges for a wide 
variety of proposal densities; however, in practice, per- 
formance is greatly improved upon judicious selection of 
this distribution. An important indicator of the perfor- 
mance of the sampling algorithm is the acceptance rate, 
the ratio between the number of accepted steps and the 
number of trial steps taken; the rule of thumb for opti- 
mal acceptance rate is c± 25% (Roberts ct al. 1997). If 
no previous information on the covariance of parameters 
(for example, from previous MCMC runs) is available, 
th e proposal dens i ty is ch osen according to the method 
of iLewis fc Bridli ()2002bf) . A random orthonormal basis 
in parameter space is chosen; the normalization to or- 
thonormal parameters relies on an initial guess on the 
width of the target posterior distribution in each direc- 
tion. Starting from the first direction of the basis vectors, 
a step of distance r in that direction is drawn from a pro- 
posal density consisting of the sum of a Gaussian and an 
exponential: 

q{r) = g (^) e(-'-/^)' + {1 - g)e-^; (5) 

Combining an exponential tail to a Gaussian proposal 
density is good practice when target probability distri- 
butions are not expected to be Gaussian. The factor 
g identifies the degre e of mixing; we chose t he value of 
5 = 2/3 suggested in ILewis fc Bridld (|2002b[) and found 
nearly optimal acceptance rate, as discussed later in the 
text. The quantity s regulates the size of a typical trial 
step with respect to the (estimated) width of the poste- 
rior in that direction, and we choose the value s = 2.4 
as sug gested by Gclma n'et al.l (jT996) and Dunklc y et al.l 
(|2005l ) . We repeat this procedure taking subsequent trial 
steps in each direction of the basis vectors, then generate 
a new random basis, and iterate this process. 

More efficient sampling is achieved if the covariances 
between the parameters are known. When a covariance 
matrix is available, eigenvectors of this matrix are used 
to select the basis of orthonormal parameters. We in- 
troduce the possibility of using an adaptive covariance 
matrix. This feature is useful in particular for param- 
eter spaces where the shape of the target distribution 
depends strongly on the position in parameter space, so 
that previously obtained covariance matrices might not 
be a good representation of the covariance in the present 
case. When this option is selected, the sampling process 
is stopped at intervals that can be chosen by the user, a 
new covariance matrix is computed on the desired frac- 
tion of the chain (for example, the second half of the 
chain; typically one would not want to use all the sam- 
ples to exclude the burn-in region), and the new matrix 
is used as a subsequent input for the proposal density. 

One of the most critical issues in MCMC sampling 
is the assessment of convergence, namely determining 
whether the probability distribution inferred from the 
chains resembles the true one. Two key requirements 
are that the entire relevant region of parameter space 
has been explored, and that the target distribution from 
the samples has become stable {i.e., it would not change 
upon further sampling). The first point is relevant in 
particular in the presence of a multi-modal distribution; 



while in general the MCMC sampling should be able to 
correctly explore such spaces, it is wise to run multi- 
ple chains starting from different locations of parame- 
ter space, and we adopt this approach. Examples of al- 
gorithms that test convergence on single and multiple 
chains respectively are the Raftery an d Lewis (hereafter 
RL) statistics (Ra fterv fc Lewij|1992[ ) and the Gelman 
and Rubin "R" test, compa ring the variance of the mean 
within and between chains ([Gelman fc Rubinlll992f) . All 
these tests are routinely performed by the publicly avail- 
able software GetDist (Lewis fc Bridle 2002a), which we 
use to analyze the chains. Convergence tests are meant to 
be performed on an unchanging proposal density. When 
the adaptive covariance matrix option is selected, the 
user can choose to stop changing the proposal density 
either when the desired acceptance rate is achieved (the 
default value is 25%), or after a certain number of consec- 
utive attempts to update the proposal density without a 
significant improvement of the acceptance rate. GalMC 
records the numbers of steps taken up to this time, so 
that they can be discarded in the computation of the 
posterior probability and in convergence tests. 

For the analysis of data on the LAEs presented in this 
work, we use the information from the RL statistics to 
discard the burn-in of the chain and to obtain indepen- 
dent samples by thinning the chain by an appropriate 
factor. To ensure that convergence has been reached, we 
require the rather stringent constraint from the R test 
R — 1 < 0.02 and follow t he other general guidelines given 
in IDunklev et all Emh . 

2.2. From parameters to observables 

To compute the likelihood associated to a point in pa- 
rameter space (in Bayesian terms, the likelihood of the 
data, given a model), we need to predict how the ob- 
served data would change as a function of the galaxies' 
physical properties that we want to measure. The param- 
eter space we have so far explored is defined by the age 
since the onset of star formation. Age, the total stellar 
mass processed into stars, Mtot, the dust reddening, de- 
fined by the excess color E(B-V), the metallicity in units 
of the Solar one, Z/Zq, the redshift, z, and the e-folding 
time T for exponential star formation history models. In 
the analysis presented in this paper the redshift has been 
fixed, since it is well determined for narrow-band selected 
LAEs. Any combination of SED parameters can be ex- 
plored, with the exception of r and metallicity. Varying 
each of these parameters requires the use of a library of 
templates, and combining them requires too high mem- 
ory usage to be convenient for an ordinary laptop com- 
puter. This does not affect our analysis significantly since 
our data do not show a strong preference for ESF vs CSF. 
Memory usage optimization for this purpose will be ex- 
plored in a subsequent paper. 

2.2.1. Stellar population synttiesis models 

As a starting point, one can use stellar population 
synthesis (SPS) models that predict the rest-frame flux 
as a function of wavelength for different ages, masses, 
metallicities, and star formation histories. We use the 
latest version of the publ icly available GALAXEV code 
(jCharlot fc Bruzuall 120101 private communication; here- 
after CBIO), although GalM C also supports the ear- 
lier version of the same code (jBruzual fc Chariot, 2003al 
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hereafter BC03). The main difference between the two 
versions is due to the different treatment of the TP-AGB 
evolution of low- and intermediate-mass stars, whose 
contribution produces redder rest-frame near-IR colors 
in the BCIO models (S. Chariot 2010, private communi- 
cation). We modified the source code of CBIO to reduce 
running time, which is essential for MCMC algorithms 
since tens of thousands of iterations are typically re- 
quired. The main changes include a significant reduction 
of output written to files by the program csp_galaxev and 
the introduction of a new analytic option for exponen- 
tially increasing star formation rate (previously available 
only through the tabulated star formation rate option, 
which is considerably slower). 

2.2.2. Star formation history 

Three classes of star formation histories (SFH) can be 
specified: simple stellar population (SSP), correspond- 
ing to an instantaneous burst; constant star formation 
rate (CSF), and exponential star formation rate (ESF) 
tp{t) — l/lrje^*/'^. In the latter case, the e-folding time 
is also a parameter of the SED fitting, and is allowed to 
be either positive or negative to represent exponentially 
decreasing or increasing star formation histories. For 
the galaxy samples considered in this paper, the allowed 
choices of star formation histories allow us to find good 
fits to the data, but in general it would be ideal to mea- 
sure the star formation history of a galaxy, rather than 
assuming a functional form for it. This method is em- 
ploye d, e.g., , by the algorithms MOPED (jHeavens et al.l 
I2OOOI ) and VESPA ()Toieiro et al.|[200l . where the star 
formation rate is recovered from the data in several age 
bins, and the data are also used to constrain the num- 
ber of underlying stellar populations. We plan to explore 
this possibility in future application of GalMC, since we 
expect that high S/N data are needed in order to obtain 
meaningful constraints on the SFH. 

We minimized the number of calls to csp_galaxev, 
which is computationally expensive, by introducing an 
"adaptive library" for exponential star formation history 
models. The user can specify the location of a library 
where the model files are stored; at each step the value 
of T is compared to those of the models available in the 
library, and csp_galaxev is run only if no close enough 
model is present. For the analysis in this paper, we have 
required an accuracy in r of min(2%, 2 Myr), but the 
criterion can be selected by the user. Each time a new 
model is computed, it is automatically stored in the li- 
brary. The final product of the call to csp_galaxev is a 
file containing the model as a function of age and star 
formation history; if required, the code also outputs a 
file containing the evolution of the stellar mass of the 
galaxy and a file containing the number of Lyman pho- 
tons as a function of time. These files are used as input 
by the program galaxev_pl, which extracts the spectrum 
at the relevant age and normalizes it to the chosen mass 
value. Thanks to these modifications, the average time 
per iteration of the GALAXEV part of our MCMC code 
is ~ 0.3 seconds on a 2.66 GHz MacBook Pro, a factor of 
^ 20 faster than the publicly available GALAXEV code. 

2.2.3. Metallicity 

The library of models available through GALAXEV 
comprise seven metallicity values, between Z/Zq = 0.005 



and Z/Zq — 5. We allow the user to select a fixed value 
of metallicity different from any of the templates, or to 
include metallicity as one of the SED fit parameters. To 
compute the model spectrum for any value of we in- 
terpolate between the two bracketing values available in 
the CBIO or BC03 stellar libraries. We tested both hn- 
ear and logarithmic interpolation, finding that logarith- 
mic interpolation (which became our method of choice) 
is more reliable. Our MCMC SED fitting runs are tech- 
nically sound; yet, some caveats need to be mentioned. 
The dependence of an SED on metallicity is by its own 
nature complicated, since different types of stars con- 
tribute in different ways and at different epochs. There- 
fore, the precision of our SED modeling is necessarily 
limited by the paucity of empirical templates, and this 
systematic uncertainty should be folded into any metal- 
licity measurement coming from photometric data rely- 
ing on the same templates. We tested the magnitude 
of this effect by interpolating between two template val- 
ues to find the predicted spectrum corresponding to one 
of the other template values, and found discrepancies of 
order 10% — 20%, significant yet not unreasonable. 

We also explored the effect of using different priors on 
the metallicity distribution. Our choice was of a uniform 
prior in log Z/Zq , motivated by the observed di stribu- 
tion for Damped Ly-a systems ([Prochaska et al.| [2003'): 
however we also performed the MCMC analysis using a 
uniform prior in Z/Zq and found a mild dependence of 
our results on the prior used for the sample with less 
signal-to-noise. This is unsurprising since in general the 
choice of priors has more influence when the data have 
less constraining power. We further discuss this issue in 
Sec. SH 

2.2.4. Impact of nebular emission 

If desired, the contribution to the model fluxes from 
nebular emission (from both continuum and lines) can 
be included. The strength of both continuum and lines is 
assumed to be proportional to the rate of H-ionizi iig pho- 
tons per second, Qo- FoUowing lSchaerer &: Vaccal ([l998), 
we describe the flux from the continuum emission as 

fx = ^^e,0o (6) 

where c is the speed of light, as is the case B recombina- 
tion coefficient for hydrogen, 1 — is the escape fraction 
of Lyman photons, and 7 is the total continuum emission 
coefficient. After choosing a nominal electron density Ng 
and temperature Te of the emitting gas, and a nomi- 
nal helium to hydrogen ratio, the emission coeflacients of 
the free-free, free-bound and two-phot ons continuu m for 
hydrogen and helium can be found in lAUeil (flOSl an d 
iFerlandl ([1980^ . Following iSchaerer fc de BarrosI (|200^ , 
we use the values e^, = 1, iVe = lOOcm'^, = lOOOOiiT, 
and [Hc/H] = -1 throughout the paper; similar values 
were al so recently measure d for a star-forming galaxy at z 
~ 2 by iRigbv et al.l ()2011f ) . The corresponding template 
is added to the reference spectrum. We have checked 
that the evolution of the number of Lyman photons and 
the total emission computed in this way are consistent 
with those output by the publicly availab le code for com- 
putation of ionizing fluxes StarB urstQQ (jLeitherer et al.l 
Il999t iVazauez fc Leithe"re^l2005D . To add the contribu- 
tion of emission lines, we assume that the luminosity of 
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the H/? line is given by 



= 4.76 X 10""e^Q, 



(7) 



and we use empirical relative line intensities for H, He, 
C, N, O, S as a function of metallicity (D. Schaerer , 
private com munication and iSchaerer fc de BarrosI 20091 : 
data from lAnders fc: Fritze-v. AlvenslebenI 120031 and 
iStorey &: Hummed 1 1995| ). This simple treatment is un- 
likely to capture the complex radiative transfer physics 
of the Lyman-a line. Since the narrow-band technique 
allows one to compute the real Lyman-a flux (as the line 
excess with respect to the continuum), LOS subtracted 
its contribution from the data and therefore we do not 
include it in our nebular emission templates. 

2.2.5. Galactic and IGM absorption 

The spectra thereby obtained can be corrected for ab- 
sorption due to dust in the galaxy and to neutral hy- 
drogen in the intergalactic medium (IGM). The former 
requires the assumption of a dust law to connect the 
SED fit parameter E(B-V) to an optical depth curve. 
GalMC currently su pports two options: the Calzetti law 
(jCalzetti et al.l[l99^ , where the parameter Ry , related to 
the size and geometry of dust grains , can be chosen by th e 
user, and a Milky- Way type law (jCardelh et al.lll989f ). 
The former is used in th e analysis of the pre sent paper, 
with a value i?„ = 4.05 (Calze tti et"aIll2000D . Starlight 
absorption by th e IGM i s taken into account using the 
prescription from iMadaul ()1995[ ). 

2.2.6. Comparison with data points 

To obtain the spectra in the observed frame we com- 
pute the luminosity distance for the cosmology specified 
by the user through the present day matter density rela- 
tive to the critical (rim.o) and the Hubble constant (Hq); 
a fiat geometry (fitot = 1) is assumed. We use the lu- 
minosity distance to convert model luminosity to fiux; 
the model fluxes are then convolved with the filter trans- 
mission curves specified by the user for each photometric 
band. Finally, the fiux densities fi, are obtained as the 
number of photons corresponding to the fluxes of the 
model divided by the reference number of photons ob- 
tained for a flat spectrum of /i^ = 1 /iJy: 



(8) 



where T*(A) is the transmission curve in the i— th band 
between Amin and Amax- The likelihood C for each model 
can now be computed (up to a normalization factor, 
which is irrelevant because we are interested in ratios 
of likelihoods) as a function of the ■ 
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where 0* is the fiux in the i-th data point and ai is its 
photometric error. 

3. ALGORITHM TESTING 

In this section we present the extensive tests we per- 
formed to ensure that our MCMC SED fit procedure pro- 
duces reliable results. This goal can be achieved by gen- 



erating mock galaxy catalogs with known physical prop- 
erties and checking that the input parameters are cor- 
rectly recovered by means of SED fitting. It is necessary 
to show not only that the estimates of the expectation 
values of parameters are unbiased, but also that the as- 
sociated uncertainties, often quoted in terms of 68% and 
95% confidence levels, are correct. 

Mock catalogs are built in the following manner. We 
use the modified BCIO SPS code to generate the spec- 
trum of a mock galaxy in 13 bands from observed-frame 
UV to observed-frame IR and follow the procedure de- 
scribed in Sec. l2.2l to obtain the corresponding fiux den- 
sity in each band. To mimic a photometric uncertainty of 
5% in each band, we add a random Gaussian noise of this 
la amplitude to the fluxes. While these assumptions on 
the uncertainty in the fiux and the extent of the available 
photometry are realistic for our stacked fluxes of LAEs 
(jGawiser et "all t2007: Lai et al.ll20"ol iGuaita et allHoiQ ') , 
the results of this test do not depend on the particular 
numbers. 

We build and test mock catalogs for both constant and 
exponential star formation histories. For the constant 
star formation case, the running time of the MCMC SED 
fitting code is limited; a chain of composed by a few thou- 
sands steps can be obtained in approximately two hours 
on a 2.66 GHz MacBook Pro laptop computer, and we 
have found that in this simple case, when three param- 
eters describe the SED, this is enough to achieve con- 
vergence. This makes the CSF scenario a suitable test 
case to check our marginalized probability distribution 
of parameters, since this task requires running GalMC 
on a large number of Mock catalogs, as explained below. 
Mock catalogs with exponential star formation history 
are used as a means to test robustness to the presence of 
degeneracies in the SED and the effect of using different 
variables and priors in describing the star formation his- 
tory. Depending on the sign and ratio of age and r, the 
timescale associated to the exponential rise or decline of 
the SFR in the mock model, the posterior probability of 
parameters can be multi-modal or even fiat if true de- 
generacies are present. By testing the MCMC SED fit in 
such scenarios, we ensured that the input parameters are 
correctly recovered even in the presence of degeneracies 
and that we don't overestimate the degree of belief in 
our results (for example, by assuming that chains have 
converged when this is not the case). Furthermore, we 
explored different parameterizations in r and were able 
to choose tFphe one that leads to the most reliable re- 
sults. 

3.1. Constant Star Formation Models: 
Parameter and error estimation recovery 

Our reference model is a galaxy characterized by con- 
stant star formation history. Age = 180 Myr, total mass 
converted into stars Mtot = 2.95 x 10® M©, and excess 
color E(B-V) = 0.2. We build 200 mock realizations of 
this galaxy, adding a Gaussian random scatter to each 
data point of la amplitude equal to the 5% photomet- 
ric error. We run the MCMC SED fitting code on the 
reference model (without scatter in the photometry) , en- 
suring that we correctly recover the input parameters, 
and hereby obtaining a reliable covariance matrix which 
gives a nearly optimal acceptance rate (around 30%). 
We then run the MCMC SED fitting code on each of 
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the mock catalogs, using this covariance matrix as an in- 
put for the proposal density. After discarding the chains 
that have not yet converged after 10000 steps, we are left 
with 172 reliable runs that we use to test the recovery of 
uncertainties. To do so, we compute the frequency with 
which each of the "true" (input) parameters is within the 
68% and 95% confidence levels from the 1-D marginal- 
ized posterior probability distributions obtained for the 
mock catalogs. The z-dimensional marginalized posterior 
distribution of parameters is obtained by integrating the 
n-dimensional posterior distribution p in n-i dimensions; 
for example, for CSF the 1-D age marginalized posterior 
probability is given by the function: 

p(Age) = J dMtot J dEp(Age,Mt„t,E). (10) 

Since a thinned MCMC chain has the property that the 
density of points is proportional to the posterior proba- 
bility (the product of likelihood and priors), the above 
integral is easily computed as a sum over the points of the 
chain that takes into account the time spent at each lo- 
cation in parameter space. By definition, we expect that 
the true values are within the region allowed at 68 (95)% 
confidence 68 (95)% of the time. We find that this is the 
case for each parameter, within the Poisson fluctuation 
error (order of 8% effect) associated with our statistics. 
We conclude that the uncertainties we report are reliable. 
We note that since we use the well-tested GetDist soft- 
ware to compute the posterior probability distributions 
and to perform convergence tests, this check ensures that 
our sampling algorithm and calculation of likelihood are 
implemented correctly. 

A more immediate visualization of the match between 
the observational errors and the corresponding uncer- 
tainties in the parameters can be obtained by plotting the 
probability distribution of the parameters in the "true" 
model (without the scatter) , and the scatter of the best- 
fit values obtained for the 172 mock catalogs. The agree- 
ment between these two quantities is however only ex- 
pected in the case of perfect (one-to-one) correspondence 
between the likelihood associated to the data for each 
model and the posterior probability, and assuming that 
the best fit represents the truth for each model. These 
assumptions are not unreasonable for these simple mod- 
els, leading to the observed agreement shown in Fig. [1] 

3.2. Exponential Star Formation Models: Parameter 
recovery 

Models with exponentially decreasing star formation 
rate have been often used in the literature, and re- 
cently increasing exponential star fo rmation rates have 
also been considered (|Maraston et aLl[2010l : iGuaita et al.l 
l2010t iLee et al1l2010( ) . While studying the constraints on 
these models is interesting, they are often affected by de- 
generacies intrinsic in the model SEDs. For this reason, 
it is necessary to pay particular attention to the use of 
suitable sampling variables and priors. We have imple- 
mented three possible choices of parametrization in r, 
using as sampling variable r, 1/|t| ln(|T|), and 1/r; all 
of these are currently available options in the code. We 
ran GalMC for each of these choices on a range of Mock 
catalogs, using different signs and numerical values for 
the ratio of age and r. The parametrization l/r was 



our final variable of choice, since we did not observe any 
bias or erroneous convergence problem for this case. This 
description has the attractive feature that the constant 
star formation rate case, which is the limit for t — > cxd for 
both positive and negative values of r, occupies a single 
spot {1/t — 0) in parameter space (unlike the other pa- 
rameterizations). This makes the interpretation of con- 
straints easier and transparent. We implemented it with 
a flat prior in Inr. In fact, during the parametrization 
selection process, we found that the use of a flat prior 
in l/r led in some cases to a poor recovery of the in- 
put parameters, as did the use of a flat prior in r when 
sampling using r as one of the parameters, which is not 
uncommon in the literature. On the other hand, no bias 
or erroneous results were found by us when using a uni- 
form prior in In |t|, as shown in Fig. [21 confirming that 
this is the best choice of prior. 

The sensitivity of the results to the choice of priors 
emphasizes the need for caution when computing con- 
straints on these models. We expect that significant im- 
provement could be achieved by using variables which 
are closer to the eigenvectors of the covariance matrix; 
we defer this treatment to a subsequent paper (Acqua- 
viva et al 2011 in prep). 

Results obtained for three illustrative mock catalogs, 
chosen to have different star formation histories (increas- 
ing and decreasing), r/Age ratios, and redshifts, are 
shown in Fig. [2j The 1-D posterior distributions are ob- 
tained by integrating the posterior probability over the 
remaining N-1 dimensions. In all the Figures, the nor- 
malization of these plots is arbitrary; we show it using 
the usual CosmoMC convention, where the height of the 
peak is one. We find good agreement, in each case within 
the region allowed at 68% confidence, between the input 
parameters and those recovered by means of SED fitting. 



3.3. Comparison with grid-based techniques 

As explained in Sec. [2l the MCMC technique allows 
one to draw samples from the posterior distribution, and 
to use them to compute any meaningful quantity associ- 
ated with it, e.g., the expectation values of parameters, 
and their uncertainties (for any quantile). This is obvi- 
ously not the only way to achieve this goal. Why should 
one make use of the MCMC technique rather than sam- 
ple the posterior distribution by means, for example, of a 
fine grid in the parameter space? There are two main ar- 
guments in favor of this choice: efficiency and reliability. 
The first can be understood as follows. In general, the 
allowed range of values for each parameter is much larger 
than the interesting range for that parameter, namely the 
region where the likelihood is significantly different from 
zero, which is the relevant one in computing integrals of 
the form given by Eq. ([l}. For example, to measure the 
age of a galaxy, one would typically want to test values 
from the youngest age that be can resolved (possibly a 
few million years) to the age of the Universe at the red- 
shift of the galaxy. Even using a logarithmic spacing, this 
spans several e-foldings. On the other hand, the width 

For Mock catalogs 1 and 3 convergence is slow and the strin- 
gent criterion R — 1 < 0.02 would not be satisfied yet. This is 
encouraging since it indicates that our criteria are indeed fairly 
conservative. 
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of the interesting region (say, to within 99% confidence 
level from the best-fit or mean age) will typically be much 
smaller (of the order of three e-foldings for the example 
of Fig. [1]). By sampling on a grid, a large fraction of 
time will be spent sampling regions of no interest, while 
in MCMC, after a burn-in period, every step, whether or 
not is accepted, is taken in the informative region. For 
an equal amount of CPU time, this means that the fi- 
nesse of the results achieved by MCMC - the accuracy 
when computing e.g., the integral in Eq. ([T]) - is much 
higher. As an example, we considered the Mock catalog 
described in 13.11 We found that in 86% of our 200 test 
runs, 10000 steps were enough to achieve convergence. 
Let us assume, to be conservative, that 20000 steps is a 
safe number of steps for a hypothetical MCMC chain to 
converge. For this three-dimensional parameter space, 
this total number corresponds to computing the model 
SED for ^ 27 values for each parameter. We performed 
the grid sampling assuming an allowed range which is 
comparable for log(Age) and E(B-V) to the top-hat pri- 
ors used for the MCMC case, and significantly narrower 
in log(Mass). The number of points falling in the relevant 
region of the posterior (shown in Fig. [T]) is only a hand- 
ful (three or four) for log(Mass) and E(B-V), and about 
ten for log(Age). Therefore, the fraction of time spent in 
the informative region is only (4 x 4 x 10)/27'^, which is 
less than 1%. This is too coarse a grid to reconstruct the 
marginalized probability distribution, and an attempt at 
it via the integration of the posterior leads to a delta 
function at the best fit found in the grid {y^ = 20, Age 
=108 Myr, Mass = 2.62 x 10^ Mq, and E(B-V) = 0.23, 
to be compared with the input values Age = 180 Myr, 
Mass = 2.95 X 10* Mq, and E(B-V) = 0.2). The situation 
only worsens for spaces with higher dimensionality, since 
the number of MCMC steps typically scales linearly in 
the number of parameters A'', while the number of steps 
in grid-based methods grows exponentially. For six pa- 
rameters, the expected numbers of MCMC steps would 
be of the order of 6 x 10^, which corresponds to sampling 
less than 10 values for each parameter on a grid in the 
same CPU time. 

Reliability is an equally important issue. The spacing 
on a grid has to be chosen arbitrarily (and, as we saw, it 
is computationally very expensive to choose a fine pace), 
and there is no guarantee that one would not miss fea- 
tures of the posterior distribution which are narrower 
than the interval between adjacent values. Conversely, 
in MCMC one can use an adaptive step size (such as 
the one implemented by us by means of the adaptive co- 
variance matrix), sampling more finely in high-likelihood 
regions. Furthermore, the use of convergence tests and 
the possibility of running multiple chains provide strong 
indications that all the interesting regions of parameter 
space have been adequately explored. These attractive 
features of MCMC technique allow for a true optimiza- 
tion of CPU time. 

4. RESULTS 

4.1. LAE samples 

As part of the MU SYC survey ()Gawiser et al.ll2006[ ). 
iGronwall et al.l (|2007D discovered a complete sample of 
162 Lyman Alpha Emitters (LAE) at 2=3.1 in a nar- 
rowband survey of the Extended Chandra Deep Field 



South (ECDF-S). The available observed- frame broad- 
band photometry for this sample encompasses six UV- 
optical bands (U, B, V, R, I, z), two near-IR bands 
(J and K), and the four IR IRAC bands. L08 elimi- 
nated from this sample 86 LAEs in regions of the IRAC 
images suffering significant contamination from nearby 
neighbors, and created samples of 18 IRAC-Detected 
and 52 IRAC-Undetected LAEs, with the "detection" 
at flux density > 0.3/i Jy equating to roughly 3(t signifi- 
cance. Six additional IRAC detections were classified as 
probable AGNs or high-dust galaxies and were analyzed 
separately. Median-stacked SEDs of these samples were 
formed and fit using BC03 models; for constant star for- 
mation rate, best-fit models showed no dust for either 
sample and stellar masses of 9 x 10^ [3 x 10*] Mq and 
ages of 1600 [160] Myr for the IRAC-Detected [Unde- 
tected] sample; the corresponding uncertainties at 68% 
confidence are reported in Table [3l The best-fit values 
obtained for exponentially declini ng SFH were within 
20% of the values obtained for CSF. lGawiser et al.l (|2007f ) 
fit a two-population SED model to the stacked IRAC- 
Undetected SED, using an exponentially declining SFH 
for each population. Although their best-fit model placed 
80% of the stellar mass in an underlying old popula- 
tion and only 20% in a 20 Myr-old starburst popu- 
lation, they were unable to rule out a single-population 
fit with age ~ 150 Myr, in good agreement with the 
results of L08. In the following, we investigate further 
the median-stacked SEDs of the L08 IRAC-Detected and 
IRAC-Undetected samples of z = 3.1 LAEs; we refer to 
these two samples as "IRAC Det" and "IRAC Und" in 
figures and tables. Our treatment differs from the pre- 
vious analyses in the use of the full probability distribu- 
tion of parameters obtained from the MCMC SED fit- 
ting code, as well as in the inclusion of the contribution 
of nebular continuum and emission lines and in the use 
of the improved CBIO SPS models. 

4.2. Physical properties and SFH 

We begin our MCMC analysis considering a four- 
dimensional parameter space defined by age, dust red- 
dening, metallicity, and stellar mass, and assuming con- 
stant star formation history. The dust reddening is 
parametrized by the excess color E(B-V) assuming a 
Calzetti dust absorption law. The input parameter of 
our modified GALAXEV code is the total mass of gas 
turned into stars Mtot (the integral of the instantaneous 
star formation rate over the age of the galaxy), but we 
report constraints on the more meaningful stellar mass 
M*, which takes into account the life cycle of stars and 
the associated mass loss. This effect is typically of order 
of 10 — 20%. We assume that none of the gas thereby 
injected into the IGM is reprocessed to stars. We use 
top-hat priors on log(Age), log(Mtot), logio(^/^o) and 
E(B-V); the allowed ranges are reported in Table [1] The 
lower bound of 10^ years for the age comes from an ed- 
ucated g uess of the applicability of the SPS models we 
consider (jBruzual &: Charlo^l2003b[ ). Our reference cos- 
mology assumes total energy density relative to critical 
^tot — 1, matter density = 0.258, and Hubble con- 
stant Hq = 73 km/sec/Mpc. In our baseline model we 
assume a Salpeter Initial Mass Function (IMF) as in L08 
and we include the contribution of nebular emission as 
described in Sec. 12.2.41 While in simple stellar popula- 
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tion (SSP) models the contribution of emission lines and 
continuum becom e negligible after ~ 20 and ~ 5 Myr 
respectively (e.g.. [Nilsson et aI]|201Cl[ ). when star forma- 
tion is ongoing Ly-a photons are continuously produced 
and this effect cannot be neglected. For example, in CSF 
models the number of Ly-a photons increases for the the 
first ~ 20 Myr and stays roughly constant thereafter. It 
is however true that for older stellar populations the rel- 
ative importance of nebular emission with respect to the 
stellar continuum decreases with time. 

Results for the IRAC Detected and Undetected sam- 
ples are presented in Fig. [3] and summarized in Table 
[21 In the table we report the mean value of parameters 
from the posterior distribution p, e.g., : 

<Age>= J dAge J dMtot j dE Age p(Age, M^^^, E) 

where the integral is over the domain specified above. 
This is easily computed by averaging over the R samples 
Vi obtained after thinning the Markov chain: 

<Age> = J dAge J dMtot J rfE Age p(Age, M^^^, E) 
^ 1 5] Age(r,) (12) 

i 

If the posterior is not symmetric (as in this case), the 
mean values can differ from the best-fit values; we re- 
turn to this issue in Sec. 14.31 We also report the 68% 
uncertainties, computed integrating the posterior from 
each side toward the high-probability region and stop- 
ping when the integral under the curve on each side 
is 16% of the total (32% in the case of E(B-V), which 
has a one-tail distribution). In Fig. [3] we show 1-D the 
marginalized posterior probability distribution for each 
parameter, and the 2-D marginalized constraints, includ- 
ing contours for the 68 and 95% confidence regions, on 
the different combinations. Degeneracies between pa- 
rameters appear here as diagonal axes of these ellipse- like 
curves; the one between age and stellar mass is evident 
in the bottom left panel. 

The different nature of the two samples is easily seen 
from Fig. [3] While both stellar populations are con- 
sistent with having no or very little dust (E(B-V) < 
0.04), the sample detected in IRAC is considerably older 
and more massive than the undetected one, with mean 
ages and stellar masses of 0.67 [0.37 - 1.18] Gyrs and 
3.2 [2.5 - 4.2] X 10^ Mq versus 6 [1-20] xlO"^ yrs and 2 
[1.1 - 3.4] xlO* M0 respectively. We are also able to set 
a constraint on the metallicity of these LAEs. For the 
IRAC detected sample, we find Z/Zq = 0.036 [0.005- 
0.07], while for the IRAC undetected sample, Z/Zq = 
0.05 [0.005-0.13]. In both cases. Solar metallicity is ex- 
cluded at more than 95% confidence. We note however 
that using a different prior (uniform in Z/Zq) slightly re- 
laxes the allowed parameter range. The metallicity con- 
straint becomes weaker for the IRAC undetected sample, 
which has lower S/N; in this case Solar metallicity is only 
excluded at 68% confidence. This interesting result is in 
alignment with previous cla ims that z=3.1 LAEs have 
metallicity lower than Solar ([Finkelstein et alll2010[ ) and 
migh t be g alaxies in their first star formation episode 
(e.q.. lHuet al..,1998. ). 



We also consider the effect of assuming an exponential 
star formation history. In this case there is an additional 
parameter of the SED fit, the e-folding time r. As ex- 
plained in Sec. 13.21 we use the variable 1 /r in the MCMC 
sampling; we apply a flat prior in log [t[ between the val- 
ues of T = -4 Gyr and r = 4 Gyr (for [t[ > 2 Gyr the SED 
is indistinguishable from that of a Constant Star Forma- 
tion History, so that the CSF case is included by this 
parametrization as the limit of l/r 0.25 Gyr^^). For 
this analysis, we fix the metallicity at the same value for 
both samples, in order to isolate the effect of assuming 
a more general SFH. We chose the value Z/Zq — 0.02, 
which is consistent with the range found earlier and for 
which we have an empirically calibrated stellar template 
available. The results are shown in Fig. |4l we consider 
the four most meaningful 2D combinations of parame- 
ters. We plot the posterior probability as function of 
1 /r to show the output of the MCMC SED fitting al- 
gorithm, but we report the constraints in terms of r for 
clarity in Table H 

In both cases the preferred value of r is close to the 
CSF value. For the Undetected sample, we find [t[ > 
0.12 Gyr at 68% confidence, while for the Detected sam- 
ple the corresponding constraint is |t| > 0.67 Gyr. This 
can perhaps be interpreted in terms of the different age 
of these two stellar populations. In fact, the SED is sen- 
sitive to the ratio Age/r; when this ratio is small, there 
aren't enough e-foldings to distinguish the effect of the 
exponential SFH from a Constant one. The change in 
the probability distribution of the other parameters is 
also shown in Fig. 2] The main effect of the inclusion 
of an extra parameter is the inability to rule out very 
old ages (leftmost panels) , so that the sharp difference in 
the age of the two populations found for CSF is some- 
what mitigated. The results for the stellar mass and dust 
reddening are stable to the change in SFH. This was ob- 
served previously by L08 for these z=3.1 samples and 
was also found to be true for redshift 2;= 2.1 LAEs by 
iGuaita et all ()2010[ ). Overall, both the IRAC Detected 
and Undetected LAEs are reasonably fit by a Constant 
Star Formation History model characterized by Age, stel- 
lar mass, and dust reddening E(B-V), with fixed metal- 
licity. The best fit models for the two cases have a of 
13.4 and and 7.6 respectively, for twelve data points and 
three free parameters. In the case of exponential SFH, 
the preferred r values are very close to the ones corre- 
sponding to CSF, but there is a modest improvement 
in the goodness-of-fit coming from adding one extra pa- 
rameter, with best-fit of 10.8 and and 5.2 respectively. 
These values are obtained for an exponentially increasing 
SFH. However, CSF is well within the range of e-folding 
times allowed at 68% confidence for both samples, and 
we conclude that the inferred physical properties of our 
LAEs are robust to different assumptions on the SFH 
and that the CSF model is a satisfying parametrization 
for both samples. 

4.3. Impact of SPS modeling 

Besides the choice of star formation history, many fea- 
tures of the stellar population synthesis modeling influ- 
ence the results of SED fitting. Our preferred model in- 
cludes nebular emission as described in Sec. 12.2.41 con- 
siders metallicity Z as one of the SED fit parameters, 
and employs the latest version of the GALAXEV code 
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and templates, CBIO. Here we change these assumptions 
one by one in a sequence and show their effect on the 
probability distribution of SED parameters and there- 
fore on the inferred properties of LAEs. We first fix the 
metallicity at the Solar value, then neglect the effect of 
nebular emission, and then use the earlier version of the 
GALAXEV code and templates, BC03. The final com- 
bination of assumptions coincides with those used in the 
analysis of LOS. 

Our results are shown in Fig. |6] and summarized in 
Table [Sj the last column also shows the best- fit , for 
twelve data points and four (for varying Z) or three free 
parameters. In all cases the inferred amount of dust red- 
dening is nearly insensitive to the assumptions on the 
SPS modeling, except for a moderate increase of the un- 
certainties when metallicity is varied, and we do not dis- 
cuss it further. This behavior is expected because we 
subtract the Ly-a fiux from our photometry before fit- 
ting. As a result, the impact of nebular emission lines 
and continuum, as well as the difference among the stellar 
templates we used, is negligible in the rest-UV region of 
the spectrum (observed frame UV-optical), whose slope 
measures dust reddening. 

The effect of assuming fixed Solar metallicity, as in 
LOS, produces an appreciable effect on the probability 
distribution of SED parameters for the IRAC Undetected 
sample. This same trend with metallicity was also ob- 
served by LOS for two discrete values of Z/Zq] however, 
a direct comparison is not possible because we use the 
CBIO models where they used the BC03 ones. The im- 
provement in the best-fit for models of low metallicity 
is especially marked for the IRAC Detected sample, as 
shown in Table [S] 

Adding the contribution of nebular emission produces 
little impact on the SED fit of the IRAC Detected sam- 
ple. In fact, this population is older, and the relative 
strength of nebular emission lines and nebular continuum 
with respect to the stellar continuum is lower. Further- 
more, the strongest emission lines (Ha, H/3 and OIII) 
mainly affect the K band; the corresponding data point 
was higher than the continuum SED, so that the SED 
parameters do not need to change to accommodate this 
feature, and instead the best-fit improves as a result. 
This behavior can be seen in the right-hand panel of Fig. 
[5l On the other hand, the inclusion of nebular emission 
has a strong effect both on the model SED and on the 
allowed parameter ranges for the younger LAEs in the 
IRAC Undetected sample, as seen in the left-hand panel 
of Fig. [5] and in Table [31 A new peak in the posterior 
probability distribution, at very young ages of ^ 1 Myr 
and very low masses around 3 xlO^ Mq, is present. As 
a result, the stellar population is interpreted as signifi- 
cantly younger and less massive when nebular emission is 
taken into account. This response of high-redshift galax- 
ies to SED fitting includi ng nebular emission is in agre e- 
ment with the findings of iSchaerer k, de BarrosI (|2009[ ). 

Finally, we consider the effect of using the BC03 stellar 
templates rather then the newer BCIO ones. Because of 
the differences in the templates discussed in Sec. 12.2.11 
when BC03 models are used, older ages and higher stellar 
masses are needed to fit the SEDs of both samples. 

The last set of assumptions in SPS modeling (no nebu- 
lar emission, solar metallicity, BC03 templates) coincides 
with the previous analysis of the same samples presented 



in LOS. This allows a direct comparison of the best-fit 
approach vs the use of mean values computed from the 
posterior distribution, shown in Fig. |6] and in Table |3l 
To compare the values of stellar masses on the same ba- 
sis, we need to account for the mass loss due to the life 
cycle of stars, which we consider and was not explic- 
itly quoted in that paper. We do so by multiplying the 
instantaneous star formation rate by the age reported 
in LOS to obtain the total mass transformed into stars 
Mtot, and transform it into stellar mass using the con- 
version factor output by our code for the same model. 
Our best-fit parameters (also shown in the figure as ma- 
genta dot-dashed vertical lines) are in exact agreement 
with those of LOS; however, the best-fit values do not lie 
exactly at the mean of the marginalized probability dis- 
tribution, leading to a moderate shift of ages and stellar 
masses. The disagreement is within the 6S% confidence 
level and is expected when the N-dimensional posterior 
distribution is asymmetric. 

5. CONCLUSIONS 

We built a Markov Chain Monte Carlo code for SED 
fitting, GalMC, designed to determine physical proper- 
ties of galaxies including age, stellar mass, dust con- 
tent, metallicity, star formation history and redshift (the 
latter was fixed in the present analysis). The purpose 
of MCMC codes is to recover the probability distribu- 
tion function (PDF) of the fit parameters by building a 
number of independent samples with the same statisti- 
cal properties as the PDF. In the Bayesian approach to 
statistical inference the relevant probability distribution 
is the posterior probability, which depends on the like- 
lihood of the parameters given the underlying data and 
on our prior beliefs about the model parameters. Sam- 
ples are built by exploring the parameter space in such a 
way that the density of sampled points is proportional to 
the probability distribution that we want to map. Once 
the samples have been obtained, they can be used to 
easily calculate the expectation values (i.e., means) of 
parameters and the associated uncertainties at any con- 
fidence level, since integrals in any dimensions become 
tractable sums. Because most of the CPU time is spent in 
high-likelihood, informative regions of parameter space, 
MCMC algorithms offer a substantial improvement in ef- 
ficiency with respect to methods based on mapping the 
posterior probability on a grid of reference parameter 
values. We found that this speed-up factor is of order ^ 
100 for a three-dimensional parameter space, and would 
rapidly increase for a larger number of parameters, since 
the expected scaling is linear in the number of parame- 
ters for MCMC, and exponential for grid-based methods. 
MCMC codes also offer warning fiags of unreliable results 
by means of convergence tests and allow one to choose 
the sampling step size adaptively in each direction. 

We conducted an extensive series of tests on Mock cat- 
alogs built from a variety of stellar populations to ensure 
that the input parameters and their uncertainties were 
correctly recovered by GalMC, and to investigate the ef- 
fect of using different priors. We then used the code 
to determine the physical properties of two samples of 
LAEs at redshift z = 3.1. Our analysis showed that, 
on average, the LAEs in the IRAC detected sample are 
older and more massive than their counterparts in the 
IRAC undetected sample, and that both populations are 
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essentially dust-free. Furthermore, LAEs at z = 3.1 have 
metallicities significantly lower than Solar, in agreement 
with recent spectroscopic studies of high-redshift LAEs. 
We investigated different assumptions on the star for- 
mation history of the LAEs, performing the fit with ei- 
ther constant or exponential (including increasing and 
decreasing) star formation rate, and found that the con- 
stant star formation model is favored by the data. Our 
results essentially confirm the findings of LOS about the 
distinct nature of the IRAC Detected and Undetected 
populations. However, important differences are found 
with respect to the previous analysis. The use of the im- 
proved SPS models from Chariot and Bruzual 2010 (as 
opposed to those from 2003) shifts the age of the older 
IRAC-detected population toward younger values, and 
therefore the mass toward lower values. Further changes 
are caused by our use of varying metallicity, by the in- 
clusion of the fiux from nebular continuum and emission 
lines, and by the use of the expectation values of pa- 
rameters computed from the posterior distribution rather 
than best-fit values. The quoted uncertainties also differ, 
since we use the Bayesian approach and compute them as 
the 68% quantile of the marginalized posterior distribu- 
tion; for the i-th parameter, this means that we integrate 
the posterior distribution in all but the i— th direction. 

The algorithm development and SED analysis con- 
ducted in the present work set the foundation for a range 
of future applications. Our analysis showed that not only 
the parameter expectation values, but also their uncer- 
tainties, depend on the assumptions made in the SPS 
modeling; for example, the inclusion of nebular emission 
noticeably worsens the ability to rule out young ages for 
the IRAC Undetected sample. In the present work we 
have compared two models with differing number of pa- 
rameters (a constant and an exponential star formation 
one). For the LAEs samples we studied, the MCMC 
sampling was directed toward very large values of the 
parameter r, similar to the CSF value, even when ESF 
was used, and the improvement in the quality of the fit 
due to use of one additional parameter was modest. We 
could conclude that the data do not strongly prefer ESF 
to CSF models. However, to better quantify this state- 
ment, or in general to assess whether the data favor the 
inclusion of additi onal parameters, we plan to use model 
selection ()Jeffrevsl lT961). 

SED fitting is a technique of ever-increasing impor- 
tance in astronomy, and it is essential that the tools 
used for statistical analysis keep up with the fast pace of 
the improvement in the data. Large- volume photometric 
surveys allow us to explore new and exciting directions, 
and this must be done while avoiding biasing assump- 
tions on the shape of the probability function, and while 
maximizing the accuracy in the reported uncertainties. 
Markov Chain Monte Carlo algorithms enable efficient, 
reliable estimation of parameter expectation values and 
uncertainties, are suitable for exploring parameter spaces 
of high dimensionality, and are able to reveal degenera- 
cies among parameters. GalMC is our implementation 
of this approach, and we hope it will prove useful for a 
wide range of astrophysical applications. 
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Fig. 1. — The red line shows the posterior probability distribution for the reference model described in the text. The blue histogram 
represents the distribution of best-fit values obtained from the analysis of 172 Mock data realizations, obtained convolving the reference 
fluxes with a random Gaussian scatter of amplitude equal to the photometric error of the reference model. The agreement confirms the 
reliability of our error estimation procedure. 
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Fig. 2. — Each row shows the marginalized 1-D posterior probability distribution for the SED fitting parameters for one of three illustrative 
mock catalogs of galaxies with Exponential Star Formation Rate (increasing and decreasing). As discussed in the text, multi-modality in 
the posterior is often observed for such models. The input parameters are shown as blue vertical lines; in each case the agreement between 
input and recovered parameters is good. 
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Fig. 3. — Top row: ID marginalized posterior distributions for constant star formation history models, for the IRAC-undetected (black, 
thin) and IRAC-detected (blue, thick) stacked samples. Bottom rows: 2D marginahzcd contours defining the 68 and 95% confidence 
regions on the corresponding pairs of parameters. 
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Fig. 4. — Top row: Impact of the assumption on the star formation history (constant, green, solid; exponential, magenta, dashed) on 
the posterior probability distribution of SED fit parameters. For this analysis we held the metallicity fixed at the value Z = 0.02 Zq. The 
main appreciable difference is the inability to rule out very old ages (comparable to the age of the Universe at z = 3.1) for exponential 
SFH, for both the IRAC Undetected (top) and IRAC Detected {bottom) samples. Our constraints on l/r show that for both samples the 
CSF value 1/t = is well within the region allowed at 68% confidence. 
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Fig, 5. — Observed-frame data points (black) and best-fit SEDs for constant star formation history and different modeling assumptions 
presented in Table [S] The colored diamonds show the best fit point for CSF in the K-band, the most sensitive to nebular emission lines. 
Left: IRAC Undetected sample. Right: IRAC Detected sample. Data points with negative median stacked fluxes are not shown. The 
blue solid line shows the best-fit SED from the run with varying metallicity, corresponding to a value of Z/Zq = 0.005 for the Detected 
sample and Z/Zq = 0.006 for the Undetected sample. 
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Fig. 6. — Impact of inclusion of nebular emission, metallicity, evolution in SPS templates, and best-fit vs mean values on the posterior 
probability distribution of SED fit parameters for the IRAC Undetected {top) and IRAC Detected (bottom) samples. Constant Star 
Formation History is assumed. The thin black vertical line and the grey shaded area indicate the best-fit models and 68% confidence levels 
from Lai et al 2008. The magenta dot-dashed vertical lines and curves show the best-fit models and probability distributions of our SED 
fit for the same SPS models and metallicity. This comparison shows the slight shift resulting from the use of mean values of the posterior 
distribution as opposed to best fit values. For the Undetected stack, the introduction of nebular emission causes the appearance of a second 
mode at very young ages and very low masses, which is particularly relevant at solar metallicity, where the contribution of the emission 
lines is stronger. Results are summarized in Table [3] 
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Parameter 


Symbol 


Type 


Prior 


Range 


Age since star formation begun 


Age 


Input 


flat in log(Age) 


logio(Age/yr) e [6.0-9.33] 


Total mass processed into stars 


Mtot 


Input 


flat in log(Mtot) 


logio(Mtot/MQ) e [4.0-15.0] 


Dust reddening 


E(B-V) 


Input 


flat in E(B-V) 


E(B-V) e [0.0-1.0] 


E-folding time 


T 


Input 


flat in log(r) 


T G [-4 Gyr, 4 Gyr] 


Metallicity 


Z/Zq 


Input 


flat in log(Z/Z0) 


logio(Z/Z0) e [-2.3, 0.7] 


Stellar Mass at t = Age 




Derived 







TABLE 1 

Input parameters of the MCMC code, priors used in the SED fit analysis, and derived parameters used to report results. 

The e-folding time T is always the input variable of the SPS code, but we have chosen use 1/r AS THE MCMC SAMPLING 

VARIABLE, AS DISCUSSED IN THE TEXT. 



Sample 


SFH 


Z/Zq 


Age (Gyr) 


E(B-V) 


M* (10^ M/Mq) 


r (Gyr) 


best fit x^/d.o.f. 


IRAC Det 


CSF 


0.02 


0.67(0.38 - 1.2] 


0.038[0 - 0.047] 


31 [24 - 39] 


oo 


13.4/9 




ESF 


0.02 


0.83[0.54 - 2.1] 


0.046[0 - 0.059] 


29 [23 - 38] 


-3.0[|t] > 0.67] 


10.8/8 


IR AC Und 


CSF 


0.02 


0.05 [0.04 - 0.16] 


0.049 [0 - 0.0G4] 


1.7 [0.6- 3.0] 


oo 


7.6/9 




ESF 


0.02 


0.22 [0.001 - 2.1] 


0.047 [0 - 0.059] 


2.0 [1.2 - 3.1] 


-2.7 [|r| > 0.12] 


5.2/8 



TABLE 2 

Results of SED fitting for constant and exponential star formation history. We report mean expectation values and 68% 

CONFIDENCE REGIONS OF PARAMETERS. FOR EXPONENTIAL SFH THE PARAMETER T IS CLOSE TO THE CSF VALUE (DEFINED BY |r| = 4 GYR 
BECAUSE OF THE PRIORS USED IN ESF MODELS). FOR BOTH SAMPLES THERE IS A MODEST IMPROVEMENT IN THE BEST-FIT WHEN USING 
EXPONENTIAL SFH, BUT CSF IS INCLUDED IN THE RANGE OF VALUES ALLOWED AT 68% CONFIDENCE, SHOWING THAT THE DATA DO NOT 

FAVOR ESF SIGNIFICANTLY. 



Sample 


Neb 


SPS 


Z/Z-o 




Age (Cyr) 


E(B-V) 


M* (10» M/Mp) 


best fit x'^/d.o f. 


IRAC Det 


Y 


CBIO 


().036[0.005 


- 0.07] 


0.(i7[0.37 - 1.18] 


0.046^0 - 0.05] 


32 [25 - 42] 


13.2/8 




Y 


CBIO 


1 




0.47[0.24 - 0.94] 


0.03910 - 0.041] 


29 [21 - 42] 


18.6/9 




N 


CBIO 


1 




0.72[0.33 - 1.3] 


0.02810 - 0.035] 


36[28 - 47] 


19.6/9 




N 


BC03 


1 




1.0[0.69 - 1.6] 


0.02510 - 0.03] 


55 [43 - 71] 


15.5/9 


LOS Det 


N 


BC03 


1 




1.6 [1.2 - 2.0] 


10 - 0.1] 


64 [43 - 85] 




IRAC Und 


Y 


CBIO 


0.05 [0.005 


- 0.13] 


0.06 [0.01 - 0.2] 


0.066 [0 - 0.085] 


2.0 [1.1- 3.4] 


6.84/8 




Y 


CBIO 


1 




0.007 [0.001 - O.OIS] 


0.035 [0 - 0.045] 


0.47 [0.26 - 0.98] 


9.92/9 




N 


CBIO 


1 




0.049 [0.021 - 0.12] 


0.04 [0 - 0.056] 


1.4 [0.89- 2.4] 


11.7/9 




N 


BC03 


1 




0.052 10.02 - 0.13] 


0.039 [0 - 0.049] 


1.5 [0.92 - 2.7] 


11.8/9 


LOS Und 


N 


BC03 


1 




0.16 [0.05 - 0.27] 


10 - 0.1] 


2.6 [0.83 - 5.7] 





TABLE 3 

Mean expectation values and 68% confidence regions from SED fitting for constant star formation history, using 

DIFFERENT ASSUMPTIONS ON INCLUSION OF NEBULAR EMISSION, SPS MODELING AND METALLICITY. SOLAR METALLICITY IS EXCLUDED AT 
95% CONFIDENCE BY THE DATA; SEE TEXT FOR MORE DETAILS. FOR COMPARISON, WE ALSO REPORT BEST-FIT VALUES AND 68% CONFIDENCE 

REGIONS OBTAINED FOR THE SAME SAMPLES BY LOS. 



